// black.h - Black option forward pricing model.
#pragma once
#define _USE_MATH_DEFINES
#include <cmath>
#include "xll_roots.h"

namespace prob {

	inline double normal_cdf(double x)
	{
		return (1 + erf(x/M_SQRT2))/2;
	}

} // namespace prob

namespace black {

	// P(F <= k) where F = f exp(-s^2t/2 + s B_t)
	// B_t/sqrt(t) <= (log k/f - s^2t/2)/s sqrt(t)
	inline double P(double f, double s, double k, double t)
	{
		auto srt = s*sqrt(t);

		return prob::normal_cdf(log(k/f)/srt + srt/2);
	}

	// P(F exp(s^2t) <= k) where F = f exp(-s^2t/2 + s B_t)
	inline double P_(double f, double s, double k, double t)
	{
		auto srt = s*sqrt(t);

		return prob::normal_cdf(log(k/f)/srt - srt/2);
	}

	// E max(k - F, 0) where F = f exp(-s^2t/2 + s B_t)
	inline double put_value(double f, double s, double k, double t)
	{
		return k*P(f, s, k, t) - f*P_(f, s, k, t);
	}
#pragma warning(push)
#pragma warning(disable: 100 101)
	// Calculate volatility given price.
	//!!! *Use classes from gsl::root*

	struct implied_volatility_params { double f; double p; double k; double t; };

	double v_function(double s, void * p) {
		implied_volatility_params * params = (implied_volatility_params *)p;
		return params->p - put_value(params->f, s, params->k, params->t);
	}

	inline double implied_volatility(double f, double p, double k, double t)
	{
		//auto v = [f,p,k,t](double s) { return p - put_value(f, s, k, t);};

		//!!! Find values that bracket the root.
		double s_lo = -10000.0, s_hi=10000.0;

		//!!! Create a 1-d solver from xll_roots.h
		const gsl_root_fsolver_type * T = gsl_root_fsolver_brent;
		gsl_root_fsolver * s = gsl_root_fsolver_alloc(T);

		//!!! Set up the solver.
		gsl_function F;
		implied_volatility_params params;
		params.f = f; params.p = p; params.k = k; params.t = t;
		F.function = &v_function;
		F.params = (void *)&params;
		gsl_root_fsolver_set(s, &F, s_lo, s_hi);

		//!!! solve for the implied volatility
		int iter = 0, max_iter = 10000;
		int status;
		double answer;
		do
		{
			iter++;
			status = gsl_root_fsolver_iterate(s);
			answer = gsl_root_fsolver_root(s);
			s_lo = gsl_root_fsolver_x_lower(s);
			s_hi = gsl_root_fsolver_x_upper(s);
			status = gsl_root_test_interval(s_lo, s_hi,
				0, 0.000001);

			if (status == GSL_SUCCESS)
				printf("Converged:\n");
		} while (status == GSL_CONTINUE && iter < max_iter);
		gsl_root_fsolver_free(s);
		return answer;
	}
#pragma warning(pop) 
} // namespace black
